###### Appendix: Figure C1
###### Robustness to Alternative Modeling Choices
gc()
rm(list = ls())
set.seed(12345)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # Note: if you are not using R Studio this command will not work, set WD to source file location manually

source("functions.R")
source("data/working_eba_function_current.R") 
require(pacman)
pacman::p_load(plm, sandwich, clubSandwich, lmtest, 
               ggplot2, numbers, dplyr, Hmisc, formattable, 
               htmltools, webshot, ggpubr, panelView, fixest, plyr, 
               scales, DataCombine)
require(webshot)
require(stringr)
date <- paste(data.table::tstrsplit(Sys.Date(), "-")[c(2:3, 1)], collapse = "")
load("data/cleaned/vat_panel_cleaned15Nov.RData")


#################### ALTERNATIVE SPECIFICATIONS ################
numSims <- 4943
depUse <- c("acc_vert", "v2x_corr")
depNames <- c("Vertical Accountability", "Control of Corruption Index")

## Country FE, quadratic time trend

modQuad <- eba_fit(dep.vars = depUse, 
                     dep.names = depNames,
                     samp.num = numSims, data = dfMerge, 
                     fe.vars = "ccode", 
                     add.terms = c("year", "year_sq"))

## Country FE, linear time trend

modLin <- eba_fit(dep.vars = depUse, 
                   dep.names = depNames,
                   samp.num = numSims, data = dfMerge, 
                   fe.vars = "ccode", 
                   add.terms = c("year"))


## Country and region-year FE, no time trend

modRegion <- eba_fit(dep.vars = depUse, 
                     dep.names = depNames,
                     samp.num = numSims, data = dfMerge, 
                     fe.vars = c("ccode", "region^year")) 

## Country FE, quadratic regional time trend

modQuadratic <- eba_fit(dep.vars = depUse, 
                   dep.names = depNames,
                   samp.num = numSims, data = dfMerge, 
                   add.terms = c("region:year", "region:year_sq"), 
                   fe.vars = "ccode")

## Country and year FE, quadratic regional time trend

modQuadraticcyfe <- eba_fit(dep.vars = depUse, 
                        dep.names = depNames,
                        samp.num = numSims, data = dfMerge, 
                        add.terms = c("region:year", "region:year_sq"))



## Country FE, linear regional time trend
modReglin <- eba_fit(dep.vars = depUse, 
                     dep.names = depNames,
                     samp.num = numSims, data = dfMerge, 
                     add.terms = "region:year", 
                     fe.vars = "ccode")

                   
## Figures

pFull <- ggarrange(modLin$acc_vert$plot + ggtitle("Country FE \n Linear Time Trend"), 
                   modQuad$acc_vert$plot + ggtitle("Country FE \n Quadratic Time Trend"), 
                   modReglin$acc_vert$plot + ggtitle("Country FE \n Linear Reg. Time Trend"), 
                   modRegion$acc_vert$plot + ggtitle("Country & Region-Year FE"), 
                   modQuadratic$acc_vert$plot + ggtitle("Country FE \n Quad. Reg. Time Trend"), 
                   modQuadraticcyfe$acc_vert$plot + ggtitle("Country & Year FE \n Quad. Reg. Time Trend"), 
                   modLin$v2x_corr$plot + ggtitle("Country FE \n Linear Time Trend"), 
                   modQuad$v2x_corr$plot + ggtitle("Country FE \n Quadratic Time Trend"), 
                   modReglin$v2x_corr$plot + ggtitle("Country FE \n Linear Reg. Time Trend"), 
                   modRegion$v2x_corr$plot + ggtitle("Country & Region-Year FE"), 
                   modQuadratic$v2x_corr$plot + ggtitle("Country FE \n Quad. Reg. Time Trend"), 
                   modQuadraticcyfe$v2x_corr$plot + ggtitle("Country & Year FE \n Quad. Reg. Time Trend"),
                   common.legend = T, legend = "bottom", ncol = 6, nrow = 2, 
                   labels = c("A", "", "", "", "", "", "B", "", "", "", "", ""))

pdf(file = "figures/appendix/appendix_figure_c1.pdf", 
    width = 24, height = 8)
print(pFull)
dev.off()


save(list = c("modLin", "modQuadratic", 
              "modRegion", "modQuad", "modReglin", "modQuadraticcyfe"), 
     file = "results/appendix_xnat_alt_specs.RData")
